Anomalous heat conduction in a carbon nanowire: Molecular dynamics calculations 
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Heat conduction of a real quasi-one dimensional material, the finite length carbon nanowire 
(CNW), inserted into the single- walled carbon nanotube (SWNT) has been studied by the molecu- 
lar dynamical (MD) method, in which both of the longitudinal as well as transverse motions of the 
chain atoms in the SWNT have been permitted. It is found that the thermal conductivity of the 
carbon nanowire is very high at room temperature, and diverges more likely with the chain length 
logarithmically. 
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I. INTRODUCTION 

When a small temperature gradient exists between 
two boundaries of a material, it is expected that heat 
will be transported through the material, which usually 

obeys the Fourier's law of conduction (j = —nVT)^ a 

well-known fact in three-dimensional systems, where j 
is the heat amount flowing through a unit surface per 
unit time,VT is the gradient of the temperature field 
over the material, and k is defined as the thermal con- 
ductivity. However, it is not clear whether the Fourier's 
law is still correct in the lower (one or two) dimensional 
systems, which stimulated a great interest in the past 
several years. It has been shown [1-14] that heat con- 
duction exhibits anomalous behavior in some one dimen- 
sional systems. For example, in the one-dimensional (ID) 
integrable systems, such as the harmonic chain and the 
monoatomic Toda model, no temperature gradient can 
be set up. In some ID nonintegrable systems, such as 
the F-K model [2], the discrete (j)^ model [4] and the 
Lorentz model [5], the temperature gradient is uniform, 
and the heat conductivity is a constant, being indepen- 
dent of system size, which means these ID systems still 
obey Fourier's law. However, in some other anharmonic 
ID systems, like the Fermi-Pasta-Ulam (FPU) /?-model 
[6, 2], the diatomic Toda chain [7], or in ID hard-particle 
gases with alternating masses [7, 8, 9], the temperature 
gradient can be formed as ^ ~ ^ and their k scales as 
~ L", where L is the system size, and a > 0. Recently, 
Wang et al. [14] studied the anomalous thermal conduc- 
tion in ID polymer chains with a modeled Hamiltonian, 
and found three types of divergent exponent a in them, 
which are caused by different couplings between longitu- 
dinal and transverse motions. However, although a great 
deal of theoretical researches on the problem had been 
made in the past several years, there exists right now 
still a lot of controversies about the divergence behavior 
of the thermal conductance in low dimensional systems, 
and many important and fundamental questions in this 
field remain unsolved. 

But all these systems said above seem to be far from 



real physical materials. What will happen for the ther- 
mal conduction in a real low dimensional material? Does 
it follow the Fourier's law or not? All of these prob- 
lems not only stimulate fundamental research interests, 
but also have great potential applications in the thermal 
conduction of the nano-materials. When the dimensions 
of electronic devices shrink to nano-scale due to the fast 
progress in the present microelectronic technology, the 
thermal conduction problem becomes more and more im- 
portant because a significant energy should be dissipated 
in a much smaller compact space. And the divergence 
of the thermal conductance with the length in the low 
dimensional materials promises possibility of making the 
more outstanding heat dissipation nano-materials, solv- 
ing the thermal dissipation problem coming from the 
miniaturization of the electronic and optical devices. So, 
it is very interesting to study the heat conduction in a 
real ID or quasi- ID physical system. 

Recently, carbon nanotubes (CNTs) have attracted 
much attention due to their remarkable electronic, ther- 
mal and mechanical properties [15]. The diameter of 
a typical nanotube ranges from several to several tens 
angstroms, and that of the smallest one is only 3 A[16]. 
While their lengths can be several /im, and even reach 
to several mm, being much larger than their diameters. 
So, carbon nanotube can be thought as a very well ID 
system. In fact, many experiments and numerical sim- 
ulations have found that the thermal conductivity n of 
SWNTs is extremely high although there exists a dis- 
tribution of the obtained n values. For example, the ob- 
served room-temperature thermal conductivity of SWNT 
rope is about 1750^5800 W/mK [17], and for individ- 
ual multiwalled carbon nanotubes (MWNTs), this value 
is larger than 3000 W/mK [18]. Using equilibrium and 
nonequilibrium molecular dynamics simulations, recent 
numerical simulations also give similar results. Berber 
et al. [19] found that, for an isolated (10, 10) nanotube 
at room temperature, n ^ 6600 W/mK. S. Maruyama et 
al. [20] claimed that hi of (5, 5) nanotube diverges as a 
power law relation with the tube length, and got a rather 
smaller hz value of about 150 ~ 500 W/mK for the (5, 5) 
tube. Zhang and Li [21] studied three armchair SWNTs, 
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i.e., (5, 5), (10, 10) and (15, 15), and found that their k.' s 
diverge as a power law, too, with their n values of about 
700 - 2200 W/mK, higher than that in Ref.[20]. Yao et 
at. [22] also studied the same three armchair tubes, and 
found their s could diverge with their lengths, and have 
the same higher value of about 400 ~ 2500 W/mK. 

At the same time, a new type of carbon structure, 
carbon nanowires (CNWs) [23] have been discovered in 
the cathode deposits prepared by hydrogen arc discharge 
evaporation of carbon rods. The CNWs are made of ex- 
traordinarily long ID linear carbon chains consisting of 
more than 100 carbon atoms inserted into the innermost 
tube (7 Adiameter) of MWNTs. The CNW can be con- 
sidered as another good example of the real ID physical 
system. In this paper, we will investigate in detail the 
heat conduction in the CNWs and pay much our atten- 
tion to the divergence behavior of its thermal conductiv- 
ity. 

In what follows we firstly introduce the model Hamil- 
tonian and calculation method. Then in Sec. Ill we give 
the main numerical results and discuss the divergence of 
thermal conductivity in the CNWs. The paper ends with 
some concluding remarks in Sec. IV. 



II. MODEL 

As well known, a bare long carbon chain can not be 
stable, and so we suppose a chain of N carbon atoms with 
the same mass rric is inserted into a single- walled carbon 
nanotube. The interaction between chain atoms is sim- 
ulated by the Tersoff- Brenner bond order potential [24], 
and the interaction between carbon chain and outside 
nanotube is described by Lennard-Jones(LJ) potential. 



the carbon atom mass. For simplicity, we assume the 
carbon atoms on outside nanotube is distributed contin- 
uously^ which is well known as the continuum model and 
used in a lot of systems [26-37]. For example, based upon 
the same idea, L.A. Girifalco et al. developed a simple 
universal graphitic potential in their paper [37]. Now, 
following them, we take external potential as: 



U(r) 



(4) 



where r and x represent the time-dependent distances of 
the wire atom to tube axis and surface element (iS, re- 
spectively. Ua is the mean surface density of tube atoms. 
In the cylindrical coordinates, U (r) can be expressed as: 



U{t) 



x)pdOdz, 



where 



= {p cos ^ — r)^ -h sin^ -\- z'^, 



(5) 



(6) 



and p is the radius of outside tube, and z is another 
two cylinder coordinates of (iS. 

Thus, when r 7^ 0, the surface integral can be simpli- 
fied to give 



U (r) = SnpnaS 
where 
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u {x) = As 
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In our simulation, £ and a are taken as 2.41 mev and 
3.4 A[25], respectively. Then the Hamiltonian of the 
chain system can be written as 



N 



where 
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Here, In is an integral related to the hyper geometric 
function, which can be made exactly in an expanded se- 
ries, and obtained result is expressed as 



Vij = fc (rij^ aijfn (^f^^-^ + bijfA (rij^ • (3) 

Here, fc (j^ij^ is a cut-off function, /r (jij^ and 

Ia (^u) Morse type functions which represent 

the attractive and repulsive effects of the potential, re- 
spectively, and Qij and bij are two parameters describing 
bond order and bond angular effects. Full details of the 
model potential are available in the original paper of Ter- 
soff and Brenner [24]. Ui is external potential exerted by 
outside nanotube. rrii is the mass of chain atoms, i.e.. 



J-n — 
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But when 6^1, i.e., r^p, the summation in the In 
will converge very slowly. So, in that case, after taking 
some algebra technique, a more efficient expression can 
be finally obtained: 
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Although this expression is an approximate one, but 
when a is smah enough, it can give very accurate result of 
and needs only to take a few terms in its summation. 
Thus combining Eq. 7, 8, 9 with Eq. 10, we obtain an 
efficient external potential, representing in high precision 
the interaction between the chain atoms and the outside 
tube. 

In this work, we only consider armchair SWNT (5, 5) 
as the outside tube because its radius is about 3.4 A, 
which is the closest to that of the innermost tube ob- 
served experimentally [23]. And the average equilibrium 
distance between the chain atoms is set to be a = 1.84 A, 
which means there are four carbon atoms in three periods 
of outside armchair nanotube. We should mention that 
at present there are NO experimental data about the 
distance between carbon atoms in the nanowire, which is 
so selected from the commensurability between both 
periods of nanowire and outside SWNT. In fact, we can 
choose different a values to study its effect on the ther- 
mal conduction of the CNW, which is beyond scope of 
the present investigation and will be left for future study. 

We determine the heat current in a temperature gradi- 
ent by nonequilibrium molecular dynamics method. Two 
atoms at each end of the CNW are subject to heat baths 
at Tl and Th respectively, which usually can be simu- 
lated by Nose-Hoover thermostats [38] . The equations of 
motion for these four atoms are written as 



h = -CRrj^Jj. (11) 

where fi is the force applied on the ith carbon atom. 
The thermal variables ("l and (^r evolve according to the 
equations 

Q = gkBTr'. 

The number of degrees of freedom for particles in ther- 
mostats is given by ^ = 6, and r is the relaxation time 
of the heat bath. The rest of the atoms follows the equa- 
tions of motion 

f, = 7„ i = 3,... ,iV-2, (13) 
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and fixed boundary conditions are assumed for the ze- 
roth and (N-hl)th atoms (ro = (0,0,0), rjv+i = 
(0,0,(Ar+l)a)). 

We first use an eighth-order Runge-Kutta algorithm to 
solve the ordinary differential equations, which provides 
more accurate results than those of the usual fourth-order 
Runge-Kutta algorithm. The time step is chosen from 
h = 0.01 to 0.05 in the unit of 0.035267 ps. Typical total 
MD steps are taken as 10*" to 10^. And for comparison, 
we also use the velocity Verlet algorithm [39] in the same 
evolution. 

The total heat flux can be expressed as follows: 

J{t) = j^Y.^i{t)ei{t), (14) 

i 

where (t)is the time-dependent coordinate of the wire 
atom i, and Si (t) is its total energy, which contains both 
of the kinetic and the potential energies. And after in- 
troducing the force on the atom i from atom j, i.e., 

Fij = the instantaneous local heat current per par- 

drj 

tide can be expressed as: 

where rij = ri — r j is the relative distance between atom 
i and j. 

Because of the linear temperature distribution, the 
classical definition on the heat conductivity can be used, 
leading to: 

K = JN/ {Tl - Tr) (16) 

When TV ^ oc, above expression corresponds to the 
coefficient of heat conductivity of the chain under tem- 
perature T = (Tl + Tr) /2. In our MD process, Tl and 
Tr are kept as 0.03 and 0.025, which correspond to real 
290K and 348K in practice, respectively. 

The alternative way to calculate the heat conductivity 
of the chain is based on the Green-Kubo formula [40] 

«. = ^/(^(r)J(0)>d., (17) 

where J (t) is defined in Eq. 14. 

In our calculations, we found that these two definitions 
always give almost equal results (the difference between 
them never exceeds a few percent). 
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FIG. 1: The temperature profiles on the chain at Tl = 0.03 
and Tr = 0.025, with N = 64 (solid line) and 128 (dash- 
dot lines). The averages are carried over a time interval of 

10^ - 10^ The distance between CNW atoms is 1.8441. 
a). Eighth-order Runge-Kutta algorithm, b). Velocity Ver let 
algorithm. 



III. NUMERICAL RESULT AND DISCUSSIONS 

In Fig. 1 we show the temperature distribution on the 
chain, calculated by both eighth-order Runge-Kutta al- 
gorithm [Fig. 1(a)] and velocity Verlet algorithm [Fig. 
1(b)]. It is seen from Fig. 1(a) and 1(b) that the linear 
temperature distribution can be obtained with both algo- 
rithms, but there are still some small differences between 
them. The velocity Verlet algorithm is a very efficient 
one, by which, a very smooth temperature distribution 
is obtained, but in this case, the Nose-Hoover boundary 
condition on the left end of the chain with higher temper- 
ature could not be well treated unless the chain is long 
enough, which may result from the sensitivity of this al- 
gorithm to the thermostat boundary condition. So, we 
will mainly use the Runge-Kutta algorithm in this work. 



FIG. 2: Thermal conductivity of the CNW as a function of 
its length, ao — 1.844 A. a) Linear- log plot, b) Log- log plot. 
The solid lines in a) and b) are used for guiding eyes. 



and the velocity Verlet algorithm is used only as a refer- 
ence. 

The calculated thermal conductivities of CNWs with 
different lengths are shown in Fig. 2, in which two types 
of scales are shown, i.e., linear- log and log- log scale. Here 
we should explain our definition on the cross section of 
CNW. The statistical radial distribution for the motion 
of wire atoms along the direction perpendicular to the 
tube axis is calculated, and obtained result is given in 

Fig. 3, which can be fitted by / (r) = A-r • exp — (^) ^ , 

with r the radial distance to the tube axis. It is found 
that the parameter h is equal to 0.151 A, and so, the final 
cross section area for the heat transport can be expressed 
as 47r6^. 

Thus in this case, the thermal conductivity n is ob- 
tained to be 142 W/mK-1323 W/mK, for the chain 
length L from 3 nm to 188 nm, which is very high. For 
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FIG. 3: Statistical radial distribution for the motion of wire 
atoms perpendicular to the tube axis. 



FIG. 4: The length dependence of thermal conductivity of 
perfect ID carbon chain model. 



comparison, it is interesting to list the value of the 
SWNTs with their lengths in the same range, which is 
about 10^ W/mK [20], or about 10^ -10^ W/mK [21, 
22]. So, the thermal conductivity of CNW is comparable 
to that of the SWNTs, at least, not much smaller than 
that of the SWNT. For example, our calculated ther- 
mal conductivity of the CNW with 512 atoms is 1.2x10^ 
W/mK. But Zhang et al. [21] got the corresponding ther- 
mal conductivity of nearly 1.6x10^ W/mK for the (5, 5) 
SWNT having 384 layers (its length is the same as that 
of the CNW), which is on the same order as our CNW 
result. 

Now, we would like to ask a question: which type of di- 
vergence behavior does our CNW system follow? Power 
law or logarithmic law? From Fig. 2 it is clearly seen that 
when the system length is increased, the hz does NOT 
show completely a linear behavior in both of linear-log 
and log- log scales, making difficult to justify which type 
of divergent behavior, power or logarithmic law is suit- 
able to the CNW. But, comparing Fig. 2a with 2b, we 
could conclude that the hi of the CNW prefers more the 
logarithmic divergence than the power law, at least for 
the middle part of the data from about 32 to 512 atoms, 
which demonstrates the divergence behavior of the CNW 
is different from that of the SWNT, following the power 
law. Why logarithmic for the CNW? We think it is just 
the transverse motions of the carbon atoms on the CNW 
to relax its thermal conductance divergence, and lead it 
to deviate from the ID power law divergence. 

In order to see more clearly the influence of the trans- 
verse motion of the CNW, we have also studied the ther- 
mal conductivity of perfect ID carbon chains connected 
by the Tersoff-Brenner bond order potential, in which, 
their transverse motions are not permitted. The initial 
equilibrium distance between their neighboring atoms 
is set to be 1.65 A, and their cross sections are deter- 
mined as follows. As well known, the cutoff distance in 



the Tersoff-Brenner bond order potential is about 2.0 A, 
which can be approximately considered as an interaction 
range between two nearby carbon chains. So, the cross 
section of a perfect ID carbon chain can be roughly esti- 
mated to be about 3.14 A^. The length dependence of the 
thermal conductivity is shown in Fig. 4. It is seen from 
Fig. 4 that the hz diverges with chain length as, cx 
with p ^ 0.39 ±0.02, and its k is about 212 W/m K-511 
W/mK, which is comparable with the result of CNW. 

Comparison between the thermal conductivities of 
both ID carbon chain and the quasi- ID CNW clearly 
show that it is indeed the transverse motions of the car- 
bon atoms on the CNW to cause its logarithmic thermal 
conductance divergence. We should emphasize that in 
real systems, the divergence of thermal conductivity will 
not be as simple as that found in the theoretical models 
such as FPU model, which probably rests with the aspect 
ratio of the system. 

The heat flux autocorrelation function in the case of 
N=64 is also calculated and shown in Fig. 5, from which 
it is seen that after a very slow decay in about 9000 
ps, the heat flux autocorrelation function decreases to a 
small value. The similar phenomenon has been observed 
by Yao et al [22], which can be understood by the fact 
that after a long enough evolution, the final state has no 
relationship with the initial state. 

Finally, we will check the validity of ensemble average 
in this low dimensional system. First, we compare those 
evolutions starting from the different initial conditions. 
The obtained results are shown in Fig. 6. Here the low 
or high initial temperature means the initial temperature 
of every atom on the chain is set to the lower or higher 
boundary temperature, respectively. And the linear ini- 
tial temperature means the initial temperature of every 
atom is chosen based upon a linear temperature distri- 
bution between the two boundaries. All the distributions 
are calculated after a time interval ^ 2 x 10^. 
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FIG. 5: Heat flux autocorrelation function of CNW with 64 
particles. 

At the first sight, four figures in Fig. 6 seem to be 
the same. In fact there are only little differences between 
them, which means our numerical results are reasonable, 
being independent of the initial conditions. However, it 
is seen from Fig. 6 that Fig. 6(d) is the smoothest one, 
which means an average over other three can give a more 
reliable result, just like averaging over a longer time in- 
terval. Thus, we can improve our calculation efficiency 



by the following method: very different initial states are 
chosen first, from which the system evolves, and after a 
period of evolution time, an average over the different 
system evolutions starting from the different initial con- 
ditions can be made. This method can be considered as 
a high efficient parallel algorithm, by which, the highest 
acceleration coefficient can be gained. 



IV. CONCLUSION 

In this paper, the heat conduction of a finite length car- 
bon chain inserted into a (5, 5) SWNT has been studied 
by using the nonequilibrium molecular dynamics method, 
in which both longitudinal and transverse motions of 
the chain atoms are permitted. The interaction between 
chain atoms and nanotubes has been simulated by a con- 
tinuous model for the nanotube. It is found that heat 
conduction of CNW does not obey the Fourier's law, and 
its thermal conductivity logarithmicly diverges with 
CNW length L as - log (L). 
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